FigureS5N

Load Data

dataset <- read.delim("raw_data/FigureS5N.txt", stringsAsFactors = FALSE)

dataset$genotype <- gsub(" ", "", dataset$genotype )
dataset$genotype <- factor(dataset$genotype)
dataset$Experiment <- factor(rep(paste0("exp", 1:(length(dataset$genotype)/length(levels(dataset$genotype)))), each=length(unique(dataset$genotype))))

dataset$UID <- factor(paste(dataset$Experiment, dataset$genotype))

# wide format
kable(dataset, row.names = F)

genotype	NT	H2O2_5uM	H2O2_30uM	H2O2_50uM	Experiment	UID
WT	1097	575	397	106	exp1	exp1 WT
YFP-ALC1	1425	1018	740	100	exp1	exp1 YFP-ALC1
WT	1768	942	462	95	exp2	exp2 WT
YFP-ALC1	1271	931	528	150	exp2	exp2 YFP-ALC1
WT	1252	741	456	95	exp3	exp3 WT
YFP-ALC1	1354	935	624	48	exp3	exp3 YFP-ALC1

library(reshape2)
# reshape to long format
dataset <- melt(dataset, variable.name = "Treatment", value.name = "Counts")

dataset$genotype <- relevel(dataset$genotype, ref = "WT")
dataset$UID <- relevel(dataset$UID, ref = "exp1 WT")

dataset$H2O2 <- gsub("NT","1",dataset$Treatment)
dataset$H2O2 <- gsub("H2O2_|uM","",dataset$H2O2)
dataset$H2O2 <- log10(as.integer(dataset$H2O2))




dataset$Offset <- NA
for(uid in levels(dataset$UID)){
        dataset$Offset[dataset$UID == uid] <- mean(dataset$Counts[dataset$UID == uid])
}

dataset$NormCounts <- dataset$Counts / dataset$Offset



dataset$Offset2 <- NA
for(gid in levels(dataset$genotype)){
        dataset$Offset2[dataset$genotype == gid] <- mean(dataset$NormCounts[dataset$genotype == gid & dataset$H2O2 == 0])
}

dataset$NormCounts2 <- dataset$NormCounts / dataset$Offset2



# long format
kable(dataset, row.names = F)

genotype	Experiment	UID	Treatment	Counts	H2O2	Offset	NormCounts	Offset2	NormCounts2
WT	exp1	exp1 WT	NT	1097	0.000000	543.75	2.0174713	2.050234	0.9840200
YFP-ALC1	exp1	exp1 YFP-ALC1	NT	1425	0.000000	820.75	1.7362169	1.776869	0.9771216
WT	exp2	exp2 WT	NT	1768	0.000000	816.75	2.1646771	2.050234	1.0558196
YFP-ALC1	exp2	exp2 YFP-ALC1	NT	1271	0.000000	720.00	1.7652778	1.776869	0.9934767
WT	exp3	exp3 WT	NT	1252	0.000000	636.00	1.9685535	2.050234	0.9601604
YFP-ALC1	exp3	exp3 YFP-ALC1	NT	1354	0.000000	740.25	1.8291118	1.776869	1.0294017
WT	exp1	exp1 WT	H2O2_5uM	575	0.698970	543.75	1.0574713	2.050234	0.5157808
YFP-ALC1	exp1	exp1 YFP-ALC1	H2O2_5uM	1018	0.698970	820.75	1.2403290	1.776869	0.6980419
WT	exp2	exp2 WT	H2O2_5uM	942	0.698970	816.75	1.1533517	2.050234	0.5625464
YFP-ALC1	exp2	exp2 YFP-ALC1	H2O2_5uM	931	0.698970	720.00	1.2930556	1.776869	0.7277158
WT	exp3	exp3 WT	H2O2_5uM	741	0.698970	636.00	1.1650943	2.050234	0.5682739
YFP-ALC1	exp3	exp3 YFP-ALC1	H2O2_5uM	935	0.698970	740.25	1.2630868	1.776869	0.7108498
WT	exp1	exp1 WT	H2O2_30uM	397	1.477121	543.75	0.7301149	2.050234	0.3561130
YFP-ALC1	exp1	exp1 YFP-ALC1	H2O2_30uM	740	1.477121	820.75	0.9016144	1.776869	0.5074175
WT	exp2	exp2 WT	H2O2_30uM	462	1.477121	816.75	0.5656566	2.050234	0.2758985
YFP-ALC1	exp2	exp2 YFP-ALC1	H2O2_30uM	528	1.477121	720.00	0.7333333	1.776869	0.4127110
WT	exp3	exp3 WT	H2O2_30uM	456	1.477121	636.00	0.7169811	2.050234	0.3497070
YFP-ALC1	exp3	exp3 YFP-ALC1	H2O2_30uM	624	1.477121	740.25	0.8429585	1.776869	0.4744067
WT	exp1	exp1 WT	H2O2_50uM	106	1.698970	543.75	0.1949425	2.050234	0.0950831
YFP-ALC1	exp1	exp1 YFP-ALC1	H2O2_50uM	100	1.698970	820.75	0.1218398	1.776869	0.0685699
WT	exp2	exp2 WT	H2O2_50uM	95	1.698970	816.75	0.1163147	2.050234	0.0567324
YFP-ALC1	exp2	exp2 YFP-ALC1	H2O2_50uM	150	1.698970	720.00	0.2083333	1.776869	0.1172474
WT	exp3	exp3 WT	H2O2_50uM	95	1.698970	636.00	0.1493711	2.050234	0.0728556
YFP-ALC1	exp3	exp3 YFP-ALC1	H2O2_50uM	48	1.698970	740.25	0.0648430	1.776869	0.0364928

Plot Data

library(ggplot2)

# raw data
ggplot(dataset, aes(x=H2O2, y=Counts)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE, aes(colour=genotype)) +
        geom_point(aes(colour=genotype, shape=Experiment), size=2) +        
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_shape_manual(values=15:20) +
        scale_color_manual(values=c("#000000","#808000"))

# NormCounts Linear
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=genotype)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=genotype), size=2) +        
        geom_smooth(method=lm, formula = y ~ x, se=FALSE) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#808000"))

# NormCounts2 Linear
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=genotype)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=genotype), size=2) +        
        geom_smooth(method=lm, formula = y ~ x, se=FALSE) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#808000"))

# NormCounts Quadratic
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=genotype)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=genotype), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)")+
        scale_color_manual(values=c("#000000","#808000"))

# NormCounts2 Quadratic
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=genotype)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=genotype), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,2), se=FALSE) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#808000"))

# NormCounts Cubic
ggplot(dataset, aes(x=H2O2, y=NormCounts, color=genotype)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=genotype), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)")+
        scale_color_manual(values=c("#000000","#808000"))

# NormCounts2 Cubic
ggplot(dataset, aes(x=H2O2, y=NormCounts2, color=genotype)) + 
        theme_bw() +
        theme(panel.grid=element_blank(), text = element_text(size=14)) +
        geom_point(aes(colour=genotype), size=2) +        
        geom_smooth(method=lm, formula = y ~ poly(x,3), se=FALSE) +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        scale_color_manual(values=c("#000000","#808000"))

library(Cairo)

cairo_pdf("FigureS5N.pdf", width = 5, height = 4, family = "Arial")

ggplot(dataset, aes(x=H2O2, y=NormCounts2)) + 
        theme_bw() +
        theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(), 
              axis.line = element_line(colour = "black"), text = element_text(size=14),
              panel.border = element_blank(), panel.background = element_blank()) +
        geom_point(aes(colour = genotype)) +
        geom_smooth(method=lm, formula = y ~ poly(x,3), se=TRUE, aes(colour = genotype), fill='#CCCCCC') +
        #facet_grid(. ~ genotype) +
        xlab(label = "H2O2 (log10 µM)") +
        ylab(label = "Normalized Counts") +
        scale_color_manual(values=c("#000000","#808000"))

dev.off()

## quartz_off_screen 
##                 2

Models

library(MASS)
library(DHARMa)
library(lme4)
library(lmerTest)
library(bbmle)

Linear formula

fit1 <- lm(Counts ~ Experiment + H2O2*genotype, data = dataset)
print(summary(fit1))

## 
## Call:
## lm(formula = Counts ~ Experiment + H2O2 * genotype, data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -245.04 -132.06  -19.98   95.30  382.01 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           1299.866    110.649  11.748 7.10e-10 ***
## Experimentexp2          86.125     94.856   0.908    0.376    
## Experimentexp3           5.875     94.856   0.062    0.951    
## H2O2                  -686.474     81.573  -8.416 1.18e-07 ***
## genotypeYFP-ALC1        60.660    135.971   0.446    0.661    
## H2O2:genotypeYFP-ALC1   35.276    115.361   0.306    0.763    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 189.7 on 18 degrees of freedom
## Multiple R-squared:  0.8839, Adjusted R-squared:  0.8517 
## F-statistic: 27.42 on 5 and 18 DF,  p-value: 7.994e-08

cat("AIC: ", AIC(fit1))

## AIC:  326.9888

simres <- simulateResiduals(fittedModel = fit1)
plot(simres)

fit2 <- lm(NormCounts ~ H2O2*genotype, data = dataset)
print(summary(fit2))

## 
## Call:
## lm(formula = NormCounts ~ H2O2 * genotype, data = dataset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30870 -0.11268 -0.01016  0.10500  0.33774 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            1.98687    0.09414  21.105 3.86e-15 ***
## H2O2                  -1.01869    0.07987 -12.754 4.60e-11 ***
## genotypeYFP-ALC1      -0.15575    0.13314  -1.170    0.256    
## H2O2:genotypeYFP-ALC1  0.16077    0.11295   1.423    0.170    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1858 on 20 degrees of freedom
## Multiple R-squared:  0.9329, Adjusted R-squared:  0.9228 
## F-statistic: 92.68 on 3 and 20 DF,  p-value: 6.638e-12

cat("AIC: ", AIC(fit2))

## AIC:  -7.066552

simres <- simulateResiduals(fittedModel = fit2)
plot(simres)

fit3 <- lm(NormCounts2 ~ H2O2*genotype, data = dataset)
print(summary(fit3))

## 
## Call:
## lm(formula = NormCounts2 ~ H2O2 * genotype, data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.173733 -0.054959 -0.005033  0.052231  0.190077 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            0.96909    0.05038  19.234 2.27e-14 ***
## H2O2                  -0.49686    0.04275 -11.624 2.38e-10 ***
## genotypeYFP-ALC1       0.06144    0.07125   0.862    0.399    
## H2O2:genotypeYFP-ALC1  0.01404    0.06045   0.232    0.819    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09941 on 20 degrees of freedom
## Multiple R-squared:  0.9301, Adjusted R-squared:  0.9196 
## F-statistic:  88.7 on 3 and 20 DF,  p-value: 9.977e-12

cat("AIC: ", AIC(fit3))

## AIC:  -37.07284

simres <- simulateResiduals(fittedModel = fit3)
plot(simres)

fit4 <- lmer(Counts ~ H2O2*genotype + (1|UID), data = dataset)
print(summary(fit4))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ H2O2 * genotype + (1 | UID)
##    Data: dataset
## 
## REML criterion at convergence: 273.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.35434 -0.61511 -0.05378  0.58692  2.21547 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  UID      (Intercept)  3127     55.92  
##  Residual             31742    178.16  
## Number of obs: 24, groups:  UID, 6
## 
## Fixed effects:
##                       Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)            1330.53      95.89   15.97  13.875 2.51e-10 ***
## H2O2                   -686.47      76.61   16.00  -8.961 1.24e-07 ***
## genotypeYFP-ALC1         60.66     135.61   15.97   0.447    0.661    
## H2O2:genotypeYFP-ALC1    35.28     108.34   16.00   0.326    0.749    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) H2O2   gYFP-A
## H2O2        -0.774              
## gntYFP-ALC1 -0.707  0.547       
## H2O2:YFP-AL  0.547 -0.707 -0.774

cat("AIC: ", AIC(fit4))

## AIC:  285.7402

simres <- simulateResiduals(fittedModel = fit4)
plot(simres)

Quadratic formula

fit5 <- lm(Counts ~ Experiment + poly(H2O2,2)*genotype, data = dataset)
print(summary(fit5))

## 
## Call:
## lm(formula = Counts ~ Experiment + poly(H2O2, 2) * genotype, 
##     data = dataset)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -226.0 -113.0  -48.0  102.0  358.9 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       634.833     77.540   8.187  4.1e-07 ***
## Experimentexp2                     86.125     94.966   0.907    0.378    
## Experimentexp3                      5.875     94.966   0.062    0.951    
## poly(H2O2, 2)1                  -2257.814    268.605  -8.406  2.9e-07 ***
## poly(H2O2, 2)2                    121.603    268.605   0.453    0.657    
## genotypeYFP-ALC1                   94.833     77.540   1.223    0.239    
## poly(H2O2, 2)1:genotypeYFP-ALC1   116.021    379.865   0.305    0.764    
## poly(H2O2, 2)2:genotypeYFP-ALC1  -477.253    379.865  -1.256    0.227    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 189.9 on 16 degrees of freedom
## Multiple R-squared:  0.8966, Adjusted R-squared:  0.8513 
## F-statistic: 19.82 on 7 and 16 DF,  p-value: 9.13e-07

cat("AIC: ", AIC(fit5))

## AIC:  328.218

simres <- simulateResiduals(fittedModel = fit5)
plot(simres)

fit6 <- lm(NormCounts ~ poly(H2O2,2)*genotype, data = dataset)
print(summary(fit6))

## 
## Call:
## lm(formula = NormCounts ~ poly(H2O2, 2) * genotype, data = dataset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.22690 -0.11902 -0.05308  0.10590  0.31224 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1.000e+00  5.115e-02  19.552 1.42e-13 ***
## poly(H2O2, 2)1                  -3.350e+00  2.506e-01 -13.372 8.67e-11 ***
## poly(H2O2, 2)2                   1.813e-01  2.506e-01   0.724   0.4786    
## genotypeYFP-ALC1                 2.783e-16  7.233e-02   0.000   1.0000    
## poly(H2O2, 2)1:genotypeYFP-ALC1  5.288e-01  3.543e-01   1.492   0.1530    
## poly(H2O2, 2)2:genotypeYFP-ALC1 -6.474e-01  3.543e-01  -1.827   0.0843 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1772 on 18 degrees of freedom
## Multiple R-squared:  0.9451, Adjusted R-squared:  0.9298 
## F-statistic: 61.92 on 5 and 18 DF,  p-value: 1.043e-10

cat("AIC: ", AIC(fit6))

## AIC:  -7.865543

simres <- simulateResiduals(fittedModel = fit6)
plot(simres)

fit7 <- lm(NormCounts2 ~ poly(H2O2,2)*genotype, data = dataset)
print(summary(fit7))

## 
## Call:
## lm(formula = NormCounts2 ~ poly(H2O2, 2) * genotype, data = dataset)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.12770 -0.06473 -0.02589  0.05405  0.17573 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      0.48775    0.02716  17.958 6.13e-13 ***
## poly(H2O2, 2)1                  -1.63419    0.13306 -12.282 3.47e-10 ***
## poly(H2O2, 2)2                   0.08843    0.13306   0.665   0.5147    
## genotypeYFP-ALC1                 0.07504    0.03841   1.954   0.0665 .  
## poly(H2O2, 2)1:genotypeYFP-ALC1  0.04617    0.18818   0.245   0.8089    
## poly(H2O2, 2)2:genotypeYFP-ALC1 -0.35076    0.18818  -1.864   0.0787 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09409 on 18 degrees of freedom
## Multiple R-squared:  0.9436, Adjusted R-squared:  0.928 
## F-statistic: 60.28 on 5 and 18 DF,  p-value: 1.308e-10

cat("AIC: ", AIC(fit7))

## AIC:  -38.24457

simres <- simulateResiduals(fittedModel = fit7)
plot(simres)

fit8 <- lmer(Counts ~ poly(H2O2,2)*genotype + (1|UID), data = dataset)
print(summary(fit8))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(H2O2, 2) * genotype + (1 | UID)
##    Data: dataset
## 
## REML criterion at convergence: 241
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.2496 -0.5899 -0.2551  0.5390  2.0928 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  UID      (Intercept)  3255     57.05  
##  Residual             31231    176.72  
## Number of obs: 24, groups:  UID, 6
## 
## Fixed effects:
##                                 Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                       665.50      60.73     4.00  10.959 0.000394
## poly(H2O2, 2)1                  -2257.81     249.92    14.00  -9.034 3.24e-07
## poly(H2O2, 2)2                    121.60     249.92    14.00   0.487 0.634098
## genotypeYFP-ALC1                   94.83      85.88     4.00   1.104 0.331429
## poly(H2O2, 2)1:genotypeYFP-ALC1   116.02     353.45    14.00   0.328 0.747575
## poly(H2O2, 2)2:genotypeYFP-ALC1  -477.25     353.45    14.00  -1.350 0.198350
##                                    
## (Intercept)                     ***
## poly(H2O2, 2)1                  ***
## poly(H2O2, 2)2                     
## genotypeYFP-ALC1                   
## poly(H2O2, 2)1:genotypeYFP-ALC1    
## poly(H2O2, 2)2:genotypeYFP-ALC1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(H2O2,2)1 pl(H2O2,2)2 gYFP-A p(H2O2,2)1:
## pl(H2O2,2)1  0.000                                           
## pl(H2O2,2)2  0.000  0.000                                    
## gntYFP-ALC1 -0.707  0.000       0.000                        
## p(H2O2,2)1:  0.000 -0.707       0.000       0.000            
## p(H2O2,2)2:  0.000  0.000      -0.707       0.000  0.000

cat("AIC: ", AIC(fit8))

## AIC:  256.9578

simres <- simulateResiduals(fittedModel = fit8)
plot(simres)

Cubic formula

fit9 <- lm(Counts ~ Experiment + poly(H2O2,3)*genotype, data = dataset)
print(summary(fit9))

## 
## Call:
## lm(formula = Counts ~ Experiment + poly(H2O2, 3) * genotype, 
##     data = dataset)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -244.67  -65.79    5.79   39.11  340.21 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       634.833     62.149  10.215 7.18e-08 ***
## Experimentexp2                     86.125     76.117   1.131   0.2769    
## Experimentexp3                      5.875     76.117   0.077   0.9396    
## poly(H2O2, 3)1                  -2257.814    215.292 -10.487 5.17e-08 ***
## poly(H2O2, 3)2                    121.603    215.292   0.565   0.5811    
## poly(H2O2, 3)3                   -410.102    215.292  -1.905   0.0775 .  
## genotypeYFP-ALC1                   94.833     62.149   1.526   0.1493    
## poly(H2O2, 3)1:genotypeYFP-ALC1   116.021    304.468   0.381   0.7089    
## poly(H2O2, 3)2:genotypeYFP-ALC1  -477.253    304.468  -1.567   0.1393    
## poly(H2O2, 3)3:genotypeYFP-ALC1  -170.664    304.468  -0.561   0.5840    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 152.2 on 14 degrees of freedom
## Multiple R-squared:  0.9419, Adjusted R-squared:  0.9045 
## F-statistic:  25.2 on 9 and 14 DF,  p-value: 3.794e-07

cat("AIC: ", AIC(fit9))

## AIC:  318.3933

simres <- simulateResiduals(fittedModel = fit9)
plot(simres)

fit10 <- lm(NormCounts ~ poly(H2O2,3)*genotype, data = dataset)
print(summary(fit10))

## 
## Call:
## lm(formula = NormCounts ~ poly(H2O2, 3) * genotype, data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.105261 -0.038084 -0.003288  0.042566  0.114443 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      1.000e+00  2.023e-02  49.429  < 2e-16 ***
## poly(H2O2, 3)1                  -3.350e+00  9.911e-02 -33.805 2.61e-16 ***
## poly(H2O2, 3)2                   1.813e-01  9.911e-02   1.829 0.086062 .  
## poly(H2O2, 3)3                  -6.343e-01  9.911e-02  -6.400 8.79e-06 ***
## genotypeYFP-ALC1                -1.694e-16  2.861e-02   0.000 1.000000    
## poly(H2O2, 3)1:genotypeYFP-ALC1  5.288e-01  1.402e-01   3.773 0.001667 ** 
## poly(H2O2, 3)2:genotypeYFP-ALC1 -6.474e-01  1.402e-01  -4.619 0.000284 ***
## poly(H2O2, 3)3:genotypeYFP-ALC1 -1.210e-01  1.402e-01  -0.863 0.400712    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07008 on 16 degrees of freedom
## Multiple R-squared:  0.9924, Adjusted R-squared:  0.989 
## F-statistic: 296.8 on 7 and 16 DF,  p-value: 1.015e-15

cat("AIC: ", AIC(fit10))

## AIC:  -51.21038

simres <- simulateResiduals(fittedModel = fit10)
plot(simres)

fit11 <- lm(NormCounts2 ~ poly(H2O2,3)*genotype, data = dataset)
print(summary(fit11))

## 
## Call:
## lm(formula = NormCounts2 ~ poly(H2O2, 3) * genotype, data = dataset)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.052134 -0.019338 -0.001694  0.020761  0.055820 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      0.48775    0.01049  46.489  < 2e-16 ***
## poly(H2O2, 3)1                  -1.63419    0.05140 -31.794 6.87e-16 ***
## poly(H2O2, 3)2                   0.08843    0.05140   1.720 0.104632    
## poly(H2O2, 3)3                  -0.30939    0.05140  -6.019 1.79e-05 ***
## genotypeYFP-ALC1                 0.07504    0.01484   5.057 0.000117 ***
## poly(H2O2, 3)1:genotypeYFP-ALC1  0.04617    0.07269   0.635 0.534263    
## poly(H2O2, 3)2:genotypeYFP-ALC1 -0.35076    0.07269  -4.825 0.000186 ***
## poly(H2O2, 3)3:genotypeYFP-ALC1 -0.11570    0.07269  -1.592 0.131009    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03634 on 16 degrees of freedom
## Multiple R-squared:  0.9925, Adjusted R-squared:  0.9893 
## F-statistic: 303.5 on 7 and 16 DF,  p-value: 8.507e-16

cat("AIC: ", AIC(fit11))

## AIC:  -82.72848

simres <- simulateResiduals(fittedModel = fit11)
plot(simres)

fit12 <- lmer(Counts ~ poly(H2O2,3)*genotype + (1|UID), data = dataset)
print(summary(fit12))

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Counts ~ poly(H2O2, 3) * genotype + (1 | UID)
##    Data: dataset
## 
## REML criterion at convergence: 204.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.57978 -0.46992  0.05654  0.30014  2.39500 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  UID      (Intercept)  7219     84.96  
##  Residual             15375    124.00  
## Number of obs: 24, groups:  UID, 6
## 
## Fixed effects:
##                                 Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                       665.50      60.73     4.00  10.959 0.000394
## poly(H2O2, 3)1                  -2257.81     175.36    12.00 -12.876  2.2e-08
## poly(H2O2, 3)2                    121.60     175.36    12.00   0.693 0.501220
## poly(H2O2, 3)3                   -410.10     175.36    12.00  -2.339 0.037475
## genotypeYFP-ALC1                   94.83      85.88     4.00   1.104 0.331430
## poly(H2O2, 3)1:genotypeYFP-ALC1   116.02     247.99    12.00   0.468 0.648277
## poly(H2O2, 3)2:genotypeYFP-ALC1  -477.25     247.99    12.00  -1.924 0.078329
## poly(H2O2, 3)3:genotypeYFP-ALC1  -170.66     247.99    12.00  -0.688 0.504421
##                                    
## (Intercept)                     ***
## poly(H2O2, 3)1                  ***
## poly(H2O2, 3)2                     
## poly(H2O2, 3)3                  *  
## genotypeYFP-ALC1                   
## poly(H2O2, 3)1:genotypeYFP-ALC1    
## poly(H2O2, 3)2:genotypeYFP-ALC1 .  
## poly(H2O2, 3)3:genotypeYFP-ALC1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(H2O2,3)1 pl(H2O2,3)2 pl(H2O2,3)3 gYFP-A p(H2O2,3)1:
## pl(H2O2,3)1  0.000                                                       
## pl(H2O2,3)2  0.000  0.000                                                
## pl(H2O2,3)3  0.000  0.000       0.000                                    
## gntYFP-ALC1 -0.707  0.000       0.000       0.000                        
## p(H2O2,3)1:  0.000 -0.707       0.000       0.000       0.000            
## p(H2O2,3)2:  0.000  0.000      -0.707       0.000       0.000  0.000     
## p(H2O2,3)3:  0.000  0.000       0.000      -0.707       0.000  0.000     
##             p(H2O2,3)2:
## pl(H2O2,3)1            
## pl(H2O2,3)2            
## pl(H2O2,3)3            
## gntYFP-ALC1            
## p(H2O2,3)1:            
## p(H2O2,3)2:            
## p(H2O2,3)3:  0.000

cat("AIC: ", AIC(fit12))

## AIC:  224.6933

simres <- simulateResiduals(fittedModel = fit12)
plot(simres)

Compare Results

ICtab(fit1,fit2,fit3,fit4,
      fit5,fit6,fit7,fit8,
      fit9,fit10,fit11,fit12,
      base=T)

##       AIC   dAIC  df
## fit11 -82.7   0.0 9 
## fit10 -51.2  31.5 9 
## fit7  -38.2  44.5 7 
## fit3  -37.1  45.7 5 
## fit6   -7.9  74.9 7 
## fit2   -7.1  75.7 5 
## fit12 224.7 307.4 10
## fit8  257.0 339.7 8 
## fit4  285.7 368.5 6 
## fit9  318.4 401.1 11
## fit1  327.0 409.7 7 
## fit5  328.2 410.9 9

Final Result

fit <- fit11

output <- coef(summary(fit))
output <- output[grep("H2O2", rownames(output)),]


rownames(output) <- gsub("poly\\(|, [1-3]\\)","", rownames(output) )
rownames(output) <- gsub("genotype",  paste0(" ",levels(dataset$genotype)[1], " vs. "), rownames(output))
rownames(output)[!(grepl("vs", rownames(output)))] <- paste(rownames(output)[!(grepl("vs", rownames(output)))], levels(dataset$genotype)[1],  sep = " in " )


# suggested result table
kable(output, row.names = T)

	Estimate	Std. Error	t value	Pr(>|t|)
H2O21 in WT	-1.6341883	0.0513987	-31.7943497	0.0000000
H2O22 in WT	0.0884280	0.0513987	1.7204317	0.1046320
H2O23 in WT	-0.3093886	0.0513987	-6.0193863	0.0000179
H2O21: WT vs. YFP-ALC1	0.0461736	0.0726887	0.6352232	0.5342627
H2O22: WT vs. YFP-ALC1	-0.3507551	0.0726887	-4.8254386	0.0001864
H2O23: WT vs. YFP-ALC1	-0.1157002	0.0726887	-1.5917212	0.1310090

write.table(output, file = "FigureS5N_Stats.txt", quote = F, sep = "\t", row.names = T, col.names = NA)

Anova

fit11a <- lm(NormCounts2 ~ poly(H2O2,3)*genotype, data = dataset)
fit11b <- lm(NormCounts2 ~ poly(H2O2,3)+genotype, data = dataset)

# anova table
anova(fit11a, fit11b)

## Analysis of Variance Table
## 
## Model 1: NormCounts2 ~ poly(H2O2, 3) * genotype
## Model 2: NormCounts2 ~ poly(H2O2, 3) + genotype
##   Res.Df      RSS Df Sum of Sq      F   Pr(>F)   
## 1     16 0.021135                                
## 2     19 0.055772 -3 -0.034637 8.7406 0.001156 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1